% 针对 Resistive (阻抗耦合)情况

clc;clear;

% 无量纲参数
param.b = 0.7;        % E_j / V0
param.r = 0.8;        % R_j / a
param.c = 0.1;        % a² C_j / L_j
param.lambda_R = 10;  % 电感耦合强度
param.n = 0.15;     % 与压电传感器相关

b = param.b;
r = param.r;
c = param.c;
n = param.n;
lambda_R = param.lambda_R;

% 外部声信号（无量纲化）
A = 0.9; M = 0; f = 0.005;

% 初始条件（无量纲）
X0 = [0.2; 0.1; 0.2; 0.1];  % [x1; y1; x2; y2]
U0 = [20; 0; 0];               % [u, v, w]

% 初始化
h = 0.005;           % 数值积分步长（与论文一致）
total_time = 1000;   % 总时间
T = 0:h:total_time;  % 正确时间向量
n_steps = length(T);
X = zeros(n_steps, length(X0)); % 阻抗耦合系统的状态
U = zeros(n_steps, length(U0)); % 洛伦兹系统的状态
X(1, :) = X0; % 初始状态
U(1, :) = U0; % 初始状态

% 混沌电流（使用洛伦兹系统生成）
param.sigma = 10; param.rho = 28; param.beta = 8/3;
sigma = param.sigma; rho= param.rho; beta = param.beta; 

% param.iChaos = 0;
% I_Chaos = param.iChaos;

% 循环更新每个时间步
for i = 1: n_steps-1
    t = T(i); % 当前时间
    u = U(i, 1);

    param.iChaos = M * u; % 混沌电流，基于洛伦兹系统
    I_Chaos = param.iChaos;  
        
    VPC = A * cos(2 * pi * f * t);
    % 更新阻抗耦合（使用rk4）
    
    X(i + 1, :) = rk4_resistive(X(i, :)', h, b, r, c, n, lambda_R, I_Chaos, VPC);

    % 更新洛伦兹系统（使用rk4）
    U(i+1, :) = rk4_lorenz(U(i, :), h, sigma, rho, beta);
end

% 绘图：x1（无量纲膜电位）与外部声信号  
figure;
subplot(2,1,1);
plot(T, X(:,1), 'k', 'LineWidth', 1.5); 


subplot(2,1,2);
plot(T, A * cos(2 * pi * f * T), 'r--', 'LineWidth', 1);
% 
% title('(a) 电感耦合下混沌共振');
